**************************************************************************
** Svolik.do								
** Joseph Wright & Daehee Bak				
** Date created: November 19, 2015						
** Date updated: December 9, 2015
**									
** Parent file: 							
**     AutocraticInstability.do        					
**									
** Using data:								
**     SvolikLeaders.dta (Authoritarian Leader Spells, 1946 – 2008							
**     SvolikCoalitions.dta (Authoritarian Ruling Coalition Spells, 1946 – 2008		
**	 GWFglobal.dta					
**************************************************************************

*****************************************
********** Svolik Comparison ************  
*****************************************
** Compare Autocratic Leaders and Ruling Coalitions
** Reported stat I: 55% of leaderfailures entail NO RC failure
** Reported stat II: 42% of leaderfailure-years entail NO RC failure
 

* Ruling Coalitions *
use SvolikCoalitions, clear
replace rc_end=17897 if ccode==510    /* recode from Dec 30 2008 to Dec 31 2008 so right-censored */
gen yr_s = year(rc_start)
gen yr_e = year(rc_end)
gen duration = yr_e -yr_s +1
expand duration
sort rc_id
gen n = _n
egen min =min(n), by(rc_id)
gen year = yr_s  if n==min
replace year = year[_n-1] + 1 if year==. & rc_id ==rc_id[_n-1]
gen rc_fail = year == yr_e 
recode rc_fail (1=0) if  rc_end_exit=="."
gen rc_fail_type =  rc_end_exit if rc_fail ==1
drop yr_* n min
rename duration rc_duration
tab rc_fail   /* 285 */
tab rc_fail  if rc_end~=17897 /* 284 RC failures */
egen m_rc_fail = max(rc_fail) if rc_end~=17897, by(ccode year)
egen tag = tag(ccode year)  if rc_end~=17897
tab m_rc_fail if tag==1
drop tag m_rc_fail
recode  rc_end (.= 18000) /* set censored to greater than max */
gen ddate = rc_end
sort ccode year ddate
save temp, replace    /* 284 RC failures in 254 country-years */

* Autocratic Leaders *
use SvolikLeaders, clear
gen stdate =date(startdate, "MDY", 2000)
gen endate = date(enddate,  "MDY", 2000)
gen yr_s = year(stdate)
gen yr_e = year(endate)
gen duration = yr_e -yr_s +1
expand duration
sort leadid
gen n = _n
egen min =min(n), by(leadid)
gen year = yr_s  if n==min
replace year = year[_n-1] + 1 if year==. & leadid ==leadid[_n-1]
gen leader_fail = year == yr_e 
gen leader_fail_type =   exit if leader_fail ==1
drop yr_* n min
rename duration leader_dduration
tab leader_fail /* 698 */
tab leader_fail if endate~=17897
egen m_leader_fail = max(leader_fail) if endate~=17897, by(ccode year)
egen tag = tag(ccode year) if endate~=17897
tab m_leader_fail if tag==1
drop tag m_leader_fail
recode  endate (.= 18000) /* set censored to greater than max */
gen ddate =endate
sort ccode year ddate  /* 635 leader failures in 528 country-years */

* Merge RC and leader data and compare *
merge ccode year ddate using temp
tab _merge
rename _merge merge1
gen rc_leader_fail = leader_fail
sort ccode year ddate
recode rc_leader_fail (1=0) if rc_end_leader~=leader
tab rc_leader_fail leader_fail, col  /* 348 Ruling coalitions */
/* 55% of leaderfailures entail NO RC failure */
tab rc_leader_fail leader_fail if endate~=17897 & rc_end~=17897, col  /* 285 Ruling coalition failures */
egen m_rc_fail= max(rc_fail) if endate~=17897 & rc_end~=17897 , by(ccode year)
egen m_leader_fail=max(leader_fail) if endate~=17897 & rc_end~=17897, by(ccode year)
egen tag = tag(ccode year) if endate~=17897 & rc_end~=17897
/* 42% of leaderfailure years entail NO RC failure */
tab m_* if tag==1,  col  
drop tag

****************************
* Svolik -- GWF comparison *
****************************
** Compare GWF regime breakdowns and Svolik Ruling Coalition transitions
** Reported stat I: 97% of GWF regime breakdown country-years entail RC failure
** Reported stat II: 86% of RC transition country-years entail GWF regime breakdown
gen cowcode = ccode
sort cowcode year
*recode m_rc_fail (1=0) if rc_end==17897  /* account for right-censored RC failures */
merge cowcode year using GWFglobal
tab _merge
rename _merge merge2
gen rc_start_yr = year(rc_start)
gen rc_end_yr = year(rc_end)
replace m_rc_fail=1 if year==rc_start_yr | year == rc_end_yr
egen max_gwf_fail = max(gwf_fail), by(cowcode year)
egen max_rc_fail = max(m_rc_fail), by(cowcode year)
egen max_leader_fail = max(m_leader_fail), by(cowcode year)
egen tag =tag(cowcode year)
keep if tag==1
recode max_rc_fail (1=0) if rc_end==17897  /* account for right-censored RC failures */
recode max_rc_fail (0=1) if cow==438 & year==2008 /* new year's eve coup in Guinea 2008 after Conte's death; not right-censored */
recode max_rc_fail (0=1) if cow==370 & year==1994 /* start of Lukashenka's regime in 1994 is a transition in GWF and Svolik */
/*
tab max_gwf_fail max_rc_fail if rc_id~=. & gwf_regime~="NA" ,row
tab max_gwf_fail max_rc_fail if rc_id~=. & gwf_regime~="NA" ,col
*/

* 5 cases of GWF regime failure but no Svolik RC transition *
tsset cow year
gen rcfaill1fl =f.max_rc_fail==1 |  max_rc_fail==1
tab  rcfaill1fl if rc_id~=. & gwf_regime~="NA" & max_gwf_fail==1  /* 97%*/
list gwf_case year rc_start_leader rc_end_leader rc_start rc_end rcfaill1fl if rc_id~=. & gwf_regime~="NA" & max_gwf_fail==1 & rcfaill1fl==0, clean
list gwf_case year rc_start_leader rc_end_leader rc_start rc_end  if rc_id~=. & gwf_regime~="NA" & max_gwf_fail==1 & rcfaill1fl==1 & max_rc_fail==0, clean

* 24 cases of Svolik RC transition but no GWF regime failure*
tsset cow year
gen gwffaill1fl =l.max_gwf_fail==1 | f.max_gwf_fail==1 |  max_gwf_fail==1
tab gwffaill1fl if rc_id~=. & gwf_regime~="NA" & gwf_regime~="" & max_rc_fail ==1,  /* 86% */
   list gwf_case year rc_start_leader rc_end_leader rc_start rc_end if rc_id~=. & gwf_regime~="NA" & gwf_regime~="" & max_rc_fail==1 & gwffaill1fl==0, clean
   list gwf_case year rc_start_leader rc_end_leader rc_start rc_end  if rc_id~=. & gwf_regime~="NA" & max_rc_fail==1 & gwffaill1fl==1 & max_gwf_fail==0, clean


   


******************************************
**********Morrison Replication************
******************************************
use temp_Polity, clear
count
merge cowcode year using svolik_merge
tab _merge
drop if _merge==2
drop _merge
tsset cowcode year


*DV* 

gen durafail1 = poldv
replace durafail1 = 0 if gwf_fail==1   /* durable failures with no GWF regime failure */

gen durafail2 = poldv
replace durafail2 = 0 if gwf_fail==0 | (gwf_fail==. & poldv==1) /* durable failures with GWF regime failure */

/* durable failure + ruling coalition failure */
gen poldv_rcfail = .
replace poldv_rcfail = durafail2 if spell_id==.
replace poldv_rcfail = 0 if poldv==0 & spell_id~=.
replace poldv_rcfail = 1 if poldv==1 & svolik_RC==1 & spell_id~=.
replace poldv_rcfail = 0 if poldv==1 & svolik_RC==0 & spell_id~=.
/* durable failure + no ruling coalition failure */
gen poldv_norcfail = .  
replace poldv_norcfail = durafail1 if spell_id==.
replace poldv_norcfail = 0 if poldv==0 & spell_id~=.
replace poldv_norcfail = 1 if poldv==1 & svolik_RC==0 & spell_id~=.
replace poldv_norcfail = 0 if poldv==1 & svolik_RC==1 & spell_id~=.


*Replication of Model 2, Table 3
logit poldv nt gdpgrowth l.lnincome d.urbanpop elf  l.lnpopdens l.pastpolitydv polityage polityageknot1 polityageknot2 if l.polity~=-77 & polity~=-77 & polity2~=-77 & l.polity2~=-77, cluster(bankscode)
estimates store m1
lroc, nograph
tab poldv if s

**************************
*** JW additions begin ***
**************************
* code Guatemala 1974 as NO ruling coalition failure b/c no GWF fail (not in Svolik sample) *
recode poldv_rcfail (.=0) if s==1 & cow==90
recode poldv_norcfail (.=1) if s==1 & cow==90

* code Indonesia 1999 as ruling coalition failure since we code multiyear Polity failures that coincide with Svolik or GWF **
recode poldv_rcfail (.=1) if s==1 & cow==850
recode poldv_norcfail (.=0) if s==1 & cow==850
************************
*** JW additions end ***
************************

*Model 2, Table 3 (durable failure + ruling coalition failure)
xi: logit poldv_rcfail nt gdpgrowth l.lnincome d.urbanpop elf  l.lnpopdens l.pastpolitydv polityage polityageknot1 polityageknot2 if l.polity~=-77 & polity~=-77 & polity2~=-77 & l.polity2~=-77, cluster(bankscode)
estimates store m2
lroc, nograph
tab poldv poldv_rcfail if e(sample)

*Model 2, Table 3 (durable failure + no ruling coalition failure)
xi: logit poldv_norcfail nt gdpgrowth l.lnincome d.urbanpop elf  l.lnpopdens l.pastpolitydv polityage polityageknot1 polityageknot2 if l.polity~=-77 & polity~=-77 & polity2~=-77 & l.polity2~=-77, cluster(bankscode)
estimates store m3
lroc, nograph
tab poldv poldv_norcfail if e(sample)



**************************
*** JW additions begin ***
**************************
*** Twice logged non-tax revenue ***  this shows that the correlation is strongest for Svolik RC failures, less so for poldv that is not Svolik RC failure
gen lnt = ln(1+abs(nt))
gen l2nt = ln(1+abs(lnt))
replace l2nt = l2nt*-1 if nt<0
hist l2nt if s==1

xi: logit poldv_rcfail l2nt gdpgrowth l.lnincome d.urbanpop elf  l.lnpopdens l.pastpolitydv polityage polityageknot1 polityageknot2 if l.polity~=-77 & polity~=-77 & polity2~=-77 & l.polity2~=-77, cluster(bankscode)
estimates store m4
lroc, nograph
xi: logit poldv_norcfail l2nt gdpgrowth l.lnincome d.urbanpop elf  l.lnpopdens l.pastpolitydv polityage polityageknot1 polityageknot2 if l.polity~=-77 & polity~=-77 & polity2~=-77 & l.polity2~=-77, cluster(bankscode)
estimates store m5
lroc, nograph
************************
*** JW additions end ***
************************
estout m1 m2 m3 using TableG1.tex, cells(b(star  fmt(%9.3f)) se(par fmt(%9.2f))) stats(ll r2 N) style(tex) replace label starlevels(* 0.10 ** 0.05)
saveold temp_Morrison_Svolik, replace



*****************
***Simulations***
*****************

use temp_Morrison_Svolik, clear
set more off
logit poldv nt gdpgrowth l.lnincome d.urbanpop elf  l.lnpopdens l.pastpolitydv polityage polityageknot1 polityageknot2 if s, cluster(bankscode)
drawnorm b1-b11, n(10000) means(e(b)) cov(e(V))  clear 
gen PROBa=.
gen PROBalow=.
gen PROBahigh=.
gen a =.
gen nt_axis = (_n-1)/50 + 0 if _n<=118
/*We want nt to run from 0 to 2.36, ie from 5pctile to 95pctile*/ 
local a=0
while `a' <= 2.36 {
	gen x_beta  = b1*`a'+b2*3.48+b3*8.06+b4*0.4+b5*0.42+b6*3.95+b7*1.86+b8*28+b9*7819+b10*33139+b11*1
	gen prob = invlogit(x_beta)
	egen probhat=mean(prob)
      _pctile prob, p(2.5,97.5) 
      scalar problow= r(r1) 
      scalar probhigh= r(r2)
 	quietly replace a = `a' 
	quietly replace PROBa =  probhat  if nt_axis==a
	quietly replace PROBalow =  problow  if nt_axis==a
	quietly replace PROBahigh =  probhigh  if nt_axis==a
	drop x_beta prob* probhat* 
	local a = `a' + .02
}
label var nt_axis "Non-tax revenue per capita (rescaled)"
label var PROBa "Replication"
sort nt_axis
save Mrep1, replace


use temp_Morrison_Svolik, clear
set more off
logit poldv_rcfail nt gdpgrowth l.lnincome d.urbanpop elf  l.lnpopdens l.pastpolitydv polityage polityageknot1 polityageknot2 if l.polity~=-77 & polity~=-77 & polity2~=-77 & l.polity2~=-77, cluster(bankscode)
drawnorm b1-b11, n(10000) means(e(b)) cov(e(V))  clear 
gen PROBb=.
gen PROBblow=.
gen PROBbhigh=.
gen a =.
gen nt_axis = (_n-1)/50 + 0 if _n<=118
/*We want nt to run from 0 to 2.36, ie from 5pctile to 95pctile*/ 
local a=0
while `a' <= 2.36 {
	gen x_beta  = b1*`a'+b2*3.48+b3*8.06+b4*0.4+b5*0.42+b6*3.95+b7*1.86+b8*28+b9*7819+b10*33139+b11*1
	gen prob = invlogit(x_beta)
	egen probhat=mean(prob)
      _pctile prob, p(2.5,97.5) 
      scalar problow= r(r1) 
      scalar probhigh= r(r2) 
 	quietly replace a = `a' 
	quietly replace PROBb =  probhat  if nt_axis==a
	quietly replace PROBblow =  problow  if nt_axis==a
	quietly replace PROBbhigh =  probhigh  if nt_axis==a
	drop x_beta prob* probhat* 
	local a = `a' + .02
}
label var nt_axis "Non-tax revenue per capita (rescaled)"
label var PROBb "Regime transition with ruling coalition failure"
sort nt_axis
save Mrep2, replace


use temp_Morrison_Svolik, clear
set more off
logit poldv_norcfail nt gdpgrowth l.lnincome d.urbanpop elf  l.lnpopdens l.pastpolitydv polityage polityageknot1 polityageknot2 if l.polity~=-77 & polity~=-77 & polity2~=-77 & l.polity2~=-77, cluster(bankscode)
drawnorm b1-b11, n(10000) means(e(b)) cov(e(V))  clear 
gen PROBc=.
gen PROBclow=.
gen PROBchigh=.
gen a =.
gen nt_axis = (_n-1)/50 + 0 if _n<=118
/*We want nt to run from 0 to 2.36, ie from 5pctile to 95pctile*/ 
local a=0
while `a' <= 2.36 {
	gen x_beta  = b1*`a'+b2*3.48+b3*8.06+b4*0.4+b5*0.42+b6*3.95+b7*1.86+b8*28+b9*7819+b10*33139+b11*1
	gen prob = invlogit(x_beta)
	egen probhat=mean(prob)
      _pctile prob, p(2.5,97.5) 
      scalar problow= r(r1) 
      scalar probhigh= r(r2) 
 	quietly replace a = `a' 
	quietly replace PROBc =  probhat  if nt_axis==a
	quietly replace PROBclow =  problow  if nt_axis==a
	quietly replace PROBchigh =  probhigh  if nt_axis==a
	drop x_beta prob* probhat* 
	local a = `a' + .02
}
label var nt_axis "Non-tax revenue per capita (rescaled)"
label var PROBc "Regime transition with ruling coalition failure"
sort nt_axis
save Mrep3, replace


merge nt_axis using Mrep1
sort nt_axis
drop _merge
merge nt_axis using Mrep2
sort nt_axis
drop _merge
merge using temp_Morrison_Svolik
set scheme lean1
label var nt "Non-tax revenue distribution"

twoway  (hist nt if nt<2.4 & s==1 & nt>0, yaxis(2)  lcolor(gs12) bin(100) scheme(lean1)) /*
*/ (line PROBa nt_axis if nt_axis<=2.4, yaxis(1))  /*
*/ (line PROBb nt_axis if nt_axis<=2.4, yaxis(1))  /*
*/ (line PROBc nt_axis if nt_axis<=2.4, yaxis(1))  /*
*/ , yscale(range (0 0.04)) ylabel(0 (.01) 0.04) xtitle("Non-tax revenue per capita (1000's)", size(3.5)) ytitle("Pr(Failure)", size(3.5)) /*
*/ legend(pos(12) col(2) ring(1) label(1 "Non-tax revenue distribution") label(2 "Replication estimate") label(3 "Ruling Coalition Failure") label(4 "Ruling Coalition Survival") )

*graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\MorrisonRep_Svolik.pdf", as(pdf) replace






****************************************
** BDM and Smith 2010 AJPS, Model 2.2 **
****************************************

use temp_BDM, clear

sort cowcode year
count
merge cowcode year using svolik_merge
tab _merge
drop _merge

*Replication model*
stset  bdm_t, fail(fail1) id(ID)
qui streg W S age Wage threat3 Wthreat NTgdp WNTgdp lGDPpcWB WlGDPpcWB  growthWB WgrowthWB , dis(wei) anc(W ) cluster(ccode)
estimates store r1
lincom NTgdp
lincom WNTgdp + NTgdp


*DV*
/* leadership failure + ruling coalition failure */
gen fail_rcfail = .
replace fail_rcfail = fail2 if spell_id==.
replace fail_rcfail = 0 if fail1==0 & spell_id~=.
replace fail_rcfail = 1 if fail1==1 & svolik_RC==1 & spell_id~=.
replace fail_rcfail = 0 if fail1==1 & svolik_RC==0 & spell_id~=.
/* leadership failure + no ruling coalition failure */
gen fail_norcfail = .  
replace fail_norcfail = fail3 if spell_id==.
replace fail_norcfail = 0 if fail1==0 & spell_id~=.
replace fail_norcfail = 1 if fail1==1 & svolik_RC==0 & spell_id~=.
replace fail_norcfail = 0 if fail1==1 & svolik_RC==1 & spell_id~=.


sum bdm_d _d fail* if e(sample)
tab fail1 fail_rcfail if e(sample)
tab fail1 fail_norcfail if e(sample)


*Only coalition failures*
stset  bdm_t, fail(fail_rcfail) id(ID)
tab _d fail_rcfail if e(sample)
qui streg W S age Wage threat3 Wthreat NTgdp WNTgdp lGDPpcWB WlGDPpcWB  growthWB WgrowthWB , dis(wei) anc(W ) cluster(cowcode)
estimates store r2
lincom NTgdp
lincom WNTgdp  + NTgdp

*Only non-coalition failures*
stset  bdm_t, fail(fail_norcfail) id(ID)
tab _d fail_norcfail if e(sample)
qui streg W S age Wage threat3 Wthreat NTgdp WNTgdp lGDPpcWB WlGDPpcWB  growthWB WgrowthWB , dis(wei) anc(W ) cluster(cowcode)
estimates store r3
lincom NTgdp
lincom WNTgdp + NTgdp

estout  r1 r2 r3 using TableG2.tex, cells(b(star  fmt(%9.3f)) se(par fmt(%9.2f))) stats(ll r2 N) style(tex) replace label starlevels(* 0.10 ** 0.05)


***Hazards***

	stset  bdm_t, fail(fail1) id(ID)
	qui streg W S age Wage threat3 Wthreat NTgdp WNTgdp lGDPpcWB WlGDPpcWB  growthWB WgrowthWB , dis(wei) cluster(ccode)
	lincom NTgdp
	lincom WNTgdp + NTgdp
	stcurve, hazard at1(NTgdp=3, WNTgdp=0, W=0) at2(NTgdp=13, WNTgdp=0, W=0) /*
	*/  legend(pos(12) col(2) ring(1) label(1 "Non-tax= 3%") label(2 "Non-tax= 13%"))/*
	*/ yscale(range (0 0.1)) xscale(range (0 12)) xlabel(0 (2) 12)  range(0 12) 
	*graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\BDMSmithRep_Svolik1.pdf", as(pdf)                            replace

	stset  bdm_t, fail(fail_rcfail) id(ID)
	qui streg W S age Wage threat3 Wthreat NTgdp WNTgdp lGDPpcWB WlGDPpcWB  growthWB WgrowthWB , dis(wei) cluster(ccode)
	lincom NTgdp
	lincom WNTgdp + NTgdp
	stcurve, hazard at1(NTgdp=3, WNTgdp=0, W=0) at2(NTgdp=13, WNTgdp=0, W=0) /*
	*/  legend(pos(12) col(2) ring(1) label(1 "Non-tax= 3%") label(2 "Non-tax= 13%"))/*
	*/ xscale(range (0 12)) xlabel(0 (2) 12)  range(0 12) 
	*graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\BDMSmithRep_Svolik2.pdf", as(pdf)                            replace

	stset  bdm_t, fail(fail_norcfail) id(ID)
	qui streg W S age Wage threat3 Wthreat NTgdp WNTgdp lGDPpcWB WlGDPpcWB  growthWB WgrowthWB , dis(wei) cluster(ccode)
	lincom NTgdp
	lincom WNTgdp + NTgdp
	stcurve, hazard at1(NTgdp=3, WNTgdp=0, W=0) at2(NTgdp=13, WNTgdp=0, W=0) /*
	*/  legend(pos(12) col(2) ring(1) label(1 "Non-tax= 3%") label(2 "Non-tax= 13%"))/*
	*/ xscale(range (0 12)) xlabel(0 (2) 12)  range(0 12) 
	*graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\BDMSmithRep_Svolik3.pdf", as(pdf)                            replace


	
	
	

***************************************
**********Ahmed Replication************
***************************************
use temp_Ahmed, clear
count
merge cowcode year using svolik_merge
tab _merge
drop if _merge==2
drop _merge
tsset cowcode year


*DV* 
/* government turnover + ruling coalition failure */
gen turnover_rcfail = .
replace turnover_rcfail = turnover2 if spell_id==.
replace turnover_rcfail = 0 if turnover==0 & spell_id~=.
replace turnover_rcfail = 1 if turnover==1 & svolik_RC==1 & spell_id~=.
replace turnover_rcfail = 0 if turnover==1 & svolik_RC==0 & spell_id~=.
/* government turnover+ no ruling coalition failure */
gen turnover_norcfail = .  
replace turnover_norcfail = turnover1 if spell_id==.
replace turnover_norcfail = 0 if turnover==0 & spell_id~=.
replace turnover_norcfail = 1 if turnover==1 & svolik_RC==0 & spell_id~=.
replace turnover_norcfail = 0 if turnover==1 & svolik_RC==1 & spell_id~=.


save temp_Ahmed_Svolik, replace


**********************************
** Ahmed 2012, Table 3 column 2 **
**********************************
*Replication
   qui probit turnover aid_remit beta beta_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy* if sumsamp==1, cluster(govtcode)
	*sum beta if e(sample), detail
	lroc, nograph
	tab turnover if e(sample)
	lincom aid_remit + beta_AR*.048   /*5 pctile*/
	lincom aid_remit + beta_AR*.5   /*95 pctile*/
	est store a1

*DPI fail and Svolik RC fail
  qui probit turnover_rcfail aid_remit beta beta_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy* if sumsamp==1, cluster(govtcode)
	*sum beta if e(sample), detail
	lroc, nograph
	tab turnover_rcfail if e(sample)
	lincom aid_remit + beta_AR*.053   /*5 pctile*/
	lincom aid_remit + beta_AR*.5   /*95 pctile*/
	est store a2

*DPI fail but No Svolik RC fail
  qui probit turnover_norcfail aid_remit beta beta_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy* if sumsamp==1, cluster(govtcode)
	*sum beta if e(sample), detail
	lroc, nograph
	tab turnover_norcfail if e(sample)
	lincom aid_remit + beta_AR*.048   /*5 pctile*/
	lincom aid_remit + beta_AR*.5   /*95 pctile*/
	est store a3

*Replication w. correctly lagged beta
  qui probit turnover aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy* if sumsamp==1, cluster(govtcode)
	*sum beta3 if e(sample), detail
	lroc, nograph
	tab turnover if e(sample)
	lincom aid_remit + beta3_AR*.048   /*5 pctile*/
	lincom aid_remit + beta3_AR*.5   /*95 pctile*/
	est store a4

*DPI fail and Svolik RC fail w. correctly lagged beta
  qui probit turnover_rcfail aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy* if sumsamp==1, cluster(govtcode)
	*sum beta3 if e(sample), detail
	lroc, nograph
	tab turnover_rcfail if e(sample)
	lincom aid_remit + beta3_AR*.053   /*5 pctile*/
	lincom aid_remit + beta3_AR*.5   /*95 pctile*/
	est store a5
 
*DPI fail but No Svolik RC fail w. correctly lagged beta
  qui probit turnover_norcfail aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy* if sumsamp==1, cluster(govtcode)
	*sum beta3 if e(sample), detail
	lroc, nograph
	tab turnover_norcfail if e(sample)
	lincom aid_remit + beta3_AR*.048   /*5 pctile*/
	lincom aid_remit + beta3_AR*.5   /*95 pctile*/
	est store a6
 
set matsize 2000
*set emptycells drop
estout a1 a2 a3 a4 a5 a6 using TableG3.tex, cells(b(star  fmt(%9.3f)) se(par fmt(%9.2f))) stats(ll r2 N) drop(time_dum* codummy* ydummy*) style(tex) replace label starlevels(* 0.10 ** 0.05)
set matsize 500

/*
 *******************
 *** Simulations ***   
 *******************
 
 * (1) *
 use temp_Ahmed_Svolik, clear
 set more off
 probit turnover aid_remit beta beta_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum* codummy* ydummy* if sumsamp==1, cluster(govtcode)
 matrix m1 = e(b)
 scalar s1 = colsof(m1)
 scalar list s1 /*343*/
 set seed 9879789
 drawnorm b1-b343, n(10000) means(e(b)) cov(e(V))  clear 
 gen PROB1dem=.
 gen PROB1demlow=.
 gen PROB1demhigh=.
 gen PROB1aut=.
 gen PROB1autlow=.
 gen PROB1authigh=.

 gen a =.
 gen ar_axis = (_n-1)/5 + 0.5 if _n<=132
 /*We want nt to run from 0.5 to 26.5, ie from 5pctile to 95pctile*/ 
 local a=0.5
 while `a' <= 26.5 {
  	gen x_betaDem  = b1*`a'+b2*.05+b3*.05*`a'+b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b343*1
  	gen x_betaAut  = b1*`a'+b2*.5+b3*.5*`a'  +b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b343*1
 	gen probD = normal(x_betaDem)
 	gen probA = normal(x_betaAut) 
 	egen probDhat = mean(probD)
      _pctile probD, p(2.5,97.5) 
      scalar probDlow= r(r1) 
      scalar probDhigh= r(r2) 
 	egen probAhat = mean(probA)
      _pctile probA, p(2.5,97.5) 
      scalar probAlow= r(r1) 
      scalar probAhigh= r(r2) 
   	quietly replace a = `a' 
 	quietly replace PROB1dem =  probDhat  if ar_axis==a
 	quietly replace PROB1demlow =  probDlow  if ar_axis==a
 	quietly replace PROB1demhigh =  probDhigh  if ar_axis==a
 	quietly replace PROB1aut =  probAhat  if ar_axis==a
 	quietly replace PROB1autlow =  probAlow  if ar_axis==a
 	quietly replace PROB1authigh =  probAhigh  if ar_axis==a
 	drop x_beta* prob*  
 	local a = `a' + .2
 }
 label var ar_axis "Aid + Remit"
 label var PROB1dem "All failures"
 label var PROB1aut "All failures"
 sort ar_axis
 twoway (line PROB1dem ar_axis)  (line PROB1aut ar_axis) 
 saveold Arep1, replace
 
  * (2) *
  use temp_Ahmed_Svolik, clear
  set more off
  probit turnover_rcfail aid_remit beta beta_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis   time_dum*  codummy* ydummy* if sumsamp==1, cluster(govtcode)
   matrix m1 = e(b)
  scalar s1 = colsof(m1)
  scalar list s1 /*343*/
  drawnorm b1-b343, n(10000) means(e(b)) cov(e(V))  clear 
 gen PROB2dem=.
 gen PROB2demlow=.
 gen PROB2demhigh=.
 gen PROB2aut=.
 gen PROB2autlow=.
 gen PROB2authigh=.
  gen a =.
  gen ar_axis = (_n-1)/5 + 0.5 if _n<=132
  /*We want nt to run from 0.5 to 26.5, ie from 5pctile to 95pctile*/ 
  local a=0.5
  while `a' <= 26.5 {
  	gen x_betaDem  = b1*`a'+b2*.05+b3*.05*`a'+b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b343*1
  	gen x_betaAut  = b1*`a'+b2*.5+b3*.5*`a'  +b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b343*1
  	gen probD = normal(x_betaDem)
  	gen probA = normal(x_betaAut) 
 	egen probDhat = mean(probD)
      _pctile probD, p(2.5,97.5) 
      scalar probDlow= r(r1) 
      scalar probDhigh= r(r2) 
 	egen probAhat = mean(probA)
      _pctile probA, p(2.5,97.5) 
      scalar probAlow= r(r1) 
      scalar probAhigh= r(r2) 
   	quietly replace a = `a' 
 	quietly replace PROB2dem =  probDhat  if ar_axis==a
 	quietly replace PROB2demlow =  probDlow  if ar_axis==a
 	quietly replace PROB2demhigh =  probDhigh  if ar_axis==a
 	quietly replace PROB2aut =  probAhat  if ar_axis==a
 	quietly replace PROB2autlow =  probAlow  if ar_axis==a
 	quietly replace PROB2authigh =  probAhigh  if ar_axis==a
  	drop x_beta* prob*  
  	local a = `a' + .2
  }
  label var ar_axis "Aid + Remit"
  label var PROB2dem "Coalition failure"
  label var PROB2aut "Coalition failure"
  sort ar_axis
  twoway (line PROB2dem ar_axis)  (line PROB2aut ar_axis) 
 saveold Arep2, replace
 
  * (3) *
  use temp_Ahmed_Svolik, clear
  set more off
  probit turnover_norcfail aid_remit beta beta_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum*  codummy* ydummy* if sumsamp==1, cluster(govtcode)
    matrix m1 = e(b)
  scalar s1 = colsof(m1)
  scalar list s1 /*343*/
  drawnorm b1-b343, n(10000) means(e(b)) cov(e(V))  clear 
  gen PROB3dem=.
  gen PROB3demlow=.
  gen PROB3demhigh=.
  gen PROB3aut=.
  gen PROB3autlow=.
  gen PROB3authigh=.
  gen a =.
  gen ar_axis = (_n-1)/5 + 0.5 if _n<=132
  /*We want nt to run from 0.5 to 26.5, ie from 5pctile to 95pctile*/ 
  local a=0.5
  while `a' <= 26.5 {
  	gen x_betaDem  = b1*`a'+b2*.05+b3*.05*`a'+b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b343*1
  	gen x_betaAut  = b1*`a'+b2*.5+b3*.5*`a'  +b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b343*1
  	gen probD = normal(x_betaDem)
  	gen probA = normal(x_betaAut) 
 	egen probDhat = mean(probD)
      _pctile probD, p(2.5,97.5) 
      scalar probDlow= r(r1) 
      scalar probDhigh= r(r2) 
 	egen probAhat = mean(probA)
      _pctile probA, p(2.5,97.5) 
      scalar probAlow= r(r1) 
      scalar probAhigh= r(r2) 
   	quietly replace a = `a' 
 	quietly replace PROB3dem =  probDhat  if ar_axis==a
 	quietly replace PROB3demlow =  probDlow  if ar_axis==a
 	quietly replace PROB3demhigh =  probDhigh  if ar_axis==a
 	quietly replace PROB3aut =  probAhat  if ar_axis==a
 	quietly replace PROB3autlow =  probAlow  if ar_axis==a
 	quietly replace PROB3authigh =  probAhigh  if ar_axis==a
  	drop x_beta* prob*  
  	local a = `a' + .2
  }
  label var ar_axis "Aid + Remit"
  label var PROB3dem "Coalition survival"
  label var PROB3aut "Coalition survival"
  sort ar_axis
  twoway (line PROB3dem ar_axis)  (line PROB3aut ar_axis) 
  saveold Arep3, replace
 
 
  * (4) *
 use temp_Ahmed_Svolik, clear
 set more off
 probit turnover aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum*  codummy* ydummy*  if sumsamp==1, cluster(govtcode)
 drawnorm b1-b343, n(10000) means(e(b)) cov(e(V))  clear 
 gen PROB4dem=.
 gen PROB4demlow=.
 gen PROB4demhigh=.
 gen PROB4aut=.
 gen PROB4autlow=.
 gen PROB4authigh=.
 gen a =.
 gen ar_axis = (_n-1)/5 + 0.5 if _n<=132
 /*We want nt to run from 0.5 to 26.5, ie from 5pctile to 95pctile*/ 
 local a=0.5
 while `a' <= 26.5 {
  	gen x_betaDem  = b1*`a'+b2*.05+b3*.05*`a'+b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b343*1
  	gen x_betaAut  = b1*`a'+b2*.5+b3*.5*`a'  +b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b343*1
 	gen probD = normal(x_betaDem)
 	gen probA = normal(x_betaAut) 
 	egen probDhat = mean(probD)
      _pctile probD, p(2.5,97.5) 
      scalar probDlow= r(r1) 
      scalar probDhigh= r(r2) 
 	egen probAhat = mean(probA)
      _pctile probA, p(2.5,97.5) 
      scalar probAlow= r(r1) 
      scalar probAhigh= r(r2) 
   	quietly replace a = `a' 
 	quietly replace PROB4dem =  probDhat  if ar_axis==a
 	quietly replace PROB4demlow =  probDlow  if ar_axis==a
 	quietly replace PROB4demhigh =  probDhigh  if ar_axis==a
 	quietly replace PROB4aut =  probAhat  if ar_axis==a
 	quietly replace PROB4autlow =  probAlow  if ar_axis==a
 	quietly replace PROB4authigh =  probAhigh  if ar_axis==a
 	drop x_beta* prob*  
 	local a = `a' + .2
 }
 label var ar_axis "Aid + Remit"
 label var PROB4dem "All failures"
 label var PROB4aut "All failures"
 sort ar_axis
 twoway (line PROB4dem ar_axis)  (line PROB4aut ar_axis) 
 saveold Arep4, replace
 
 * (5) *
 use temp_Ahmed_Svolik, clear
 set more off
 probit turnover_rcfail aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum*  codummy* ydummy* if sumsamp==1, cluster(govtcode)
   matrix m1 = e(b)
 scalar s1 = colsof(m1)
 scalar list s1 /*343*/
 drawnorm b1-b343, n(10000) means(e(b)) cov(e(V))  clear 
 gen PROB5dem=.
 gen PROB5demlow=.
 gen PROB5demhigh=.
 gen PROB5aut=.
 gen PROB5autlow=.
 gen PROB5authigh=.
 gen a =.
 gen ar_axis = (_n-1)/5 + 0.5 if _n<=132
 /*We want nt to run from 0.5 to 26.5, ie from 5pctile to 95pctile*/ 
 local a=0.5
 while `a' <= 26.5 {
  	gen x_betaDem  = b1*`a'+b2*.05+b3*.05*`a'+b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b343*1
  	gen x_betaAut  = b1*`a'+b2*.5+b3*.5*`a'  +b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b343*1
 	gen probD = normal(x_betaDem)
 	gen probA = normal(x_betaAut) 
 	egen probDhat = mean(probD)
      _pctile probD, p(2.5,97.5) 
      scalar probDlow= r(r1) 
      scalar probDhigh= r(r2) 
 	egen probAhat = mean(probA)
      _pctile probA, p(2.5,97.5) 
      scalar probAlow= r(r1) 
      scalar probAhigh= r(r2) 
   	quietly replace a = `a' 
 	quietly replace PROB5dem =  probDhat  if ar_axis==a
 	quietly replace PROB5demlow =  probDlow  if ar_axis==a
 	quietly replace PROB5demhigh =  probDhigh  if ar_axis==a
 	quietly replace PROB5aut =  probAhat  if ar_axis==a
 	quietly replace PROB5autlow =  probAlow  if ar_axis==a
 	quietly replace PROB5authigh =  probAhigh  if ar_axis==a
 	drop x_beta* prob*  
 	local a = `a' + .2
 }
 label var ar_axis "Aid + Remit"
  label var PROB5dem "Coalition failure"
  label var PROB5aut "Coalition failure"
 sort ar_axis
 twoway (line PROB5dem ar_axis)  (line PROB5aut ar_axis) 
 saveold Arep5, replace
 
  * (6) *
 use temp_Ahmed_Svolik, clear
 set more off
 probit turnover_norcfail aid_remit beta3 beta3_AR finittrm lngdpcap ggdpcap lnpop war ldis hdis time_dum*  codummy* ydummy*  if sumsamp==1, cluster(govtcode)
 drawnorm b1-b343, n(10000) means(e(b)) cov(e(V))  clear 
 gen PROB6dem=.
 gen PROB6demlow=.
 gen PROB6demhigh=.
 gen PROB6aut=.
 gen PROB6autlow=.
 gen PROB6authigh=.
 gen a =.
 gen ar_axis = (_n-1)/5 + 0.5 if _n<=132
 /*We want nt to run from 0.5 to 26.5, ie from 5pctile to 95pctile*/ 
 local a=0.5
 while `a' <= 26.5 {
  	gen x_betaDem  = b1*`a'+b2*.05+b3*.05*`a'+b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b343*1
  	gen x_betaAut  = b1*`a'+b2*.5+b3*.5*`a'  +b4*1+b5*6.86+b6*1.34+b7*16.06+b8*0+b9*0+b10*0+b343*1
 	gen probD = normal(x_betaDem)
 	gen probA = normal(x_betaAut) 
 	egen probDhat = mean(probD)
      _pctile probD, p(2.5,97.5) 
      scalar probDlow= r(r1) 
      scalar probDhigh= r(r2) 
 	egen probAhat = mean(probA)
      _pctile probA, p(2.5,97.5) 
      scalar probAlow= r(r1) 
      scalar probAhigh= r(r2) 
   	quietly replace a = `a' 
 	quietly replace PROB6dem =  probDhat  if ar_axis==a
 	quietly replace PROB6demlow =  probDlow  if ar_axis==a
 	quietly replace PROB6demhigh =  probDhigh  if ar_axis==a
 	quietly replace PROB6aut =  probAhat  if ar_axis==a
 	quietly replace PROB6autlow =  probAlow  if ar_axis==a
 	quietly replace PROB6authigh =  probAhigh  if ar_axis==a
 	drop x_beta* prob*  
 	local a = `a' + .2
 }
 label var ar_axis "Aid + Remit"
  label var PROB6dem "Coalition survival"
  label var PROB6aut "Coalition survival"
 sort ar_axis
 twoway (line PROB6dem ar_axis)  (line PROB6aut ar_axis) 
 saveold Arep6, replace
 
 
 use Arep1, clear
 sort ar_axis
 merge ar_axis using Arep2
 sort ar_axis
 drop _merge
 merge ar_axis using Arep3
 sort ar_axis
 drop _merge
 merge ar_axis using Arep4
 sort ar_axis
 drop _merge
 merge ar_axis using Arep5
 sort ar_axis
 drop _merge
 merge ar_axis using Arep6
 sort ar_axis
 drop _merge
 merge using temp_Ahmed_Svolik
 set scheme lean1
 label var aid_remit "Aid + Remit distribution"
 sort ar_axis
 drop _merge
 
twoway  (hist aid_remit if aid_remit<=26.5 & sumsam==1 & aid_remit>=0.5, lcolor(gs12) bin(100)) /*
*/ (line PROB1aut ar_axis if ar_axis<=26.5, yaxis(1))  /*
*/ (line PROB2aut ar_axis if ar_axis<=26.5, yaxis(1))  /*
*/ (line PROB3aut ar_axis if ar_axis<=26.5, yaxis(1))  /*
*/ ,  xtitle("Aid + Remit (%GDP)", size(3.5)) ytitle("Pr(Failure)", size(3.5)) /*
*/  legend(pos(12) col(2) ring(1) label(1 "Aid + Remit distribution")  label(2 "All failures") label(3 "Coalition failure") label(4 "Coalition survival"))  xscale(range (0 26)) xlabel(0 (5) 25)
*graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\AhmedReplication_Svolik.pdf", as(pdf)                            replace

twoway  (hist aid_remit if aid_remit<=26.5 & sumsam==1 & aid_remit>=0.5, lcolor(gs12) bin(100)) /*
*/ (line PROB4aut ar_axis if ar_axis<=26.5, yaxis(1))  /*
*/ (line PROB5aut ar_axis if ar_axis<=26.5, yaxis(1))  /*
*/ (line PROB6aut ar_axis if ar_axis<=26.5, yaxis(1))  /*
*/ ,  xtitle("Aid + Remit (%GDP)", size(3.5)) ytitle("Pr(Failure)", size(3.5)) /*
*/  legend(pos(12) col(2) ring(1) label(1 "Aid + Remit distribution")  label(2 "All failures") label(3 "Coalition failure") label(4 "Coalition survival"))  xscale(range (0 26)) xlabel(0 (5) 25)
*graph export "C:\Users\jwright\Documents\My Dropbox\Research\Autocratic Instability\AhmedCorrected_Svolik.pdf", as(pdf)                            replace

*/
	